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Abstract 

We present a method for simulating fluid vesicles with in-plane orientational ordering. The method involves 
computation of local curvature tensor and parallel transport of the orientational field on a randomly triangulated 
surface. It is shown that the model reproduces the known equilibrium conformation of fluid membranes and work 
well for a large range of bending rigidities. Introduction of nematic ordering leads to stiffening of the membrane. 
Nematic ordering can also result in anisotropic rigidity on the surface leading to formation of membrane tubes. 

PACS numbers: PACS-87.16.D-, Membranes, bilayers and vesicles. PACS-05.40.-a Fluctuation phenomena, random pro- 
cesses, noise and Brownian motion. PACS-05.70.Np Interfaces and surface thermodynamics 



I. INTRODUCTION 

The phenomenological models of fluid mem- 
brane conformations have a remarkable simplicity 
due to the symmetry constraints they must obey 
[1 . However, elementary questions on the large 
scale properties of fluid membranes remain unre- 
solved due to the technical complexity in analyz- 
ing the statistical mechanics of these membrane 
models. This is in particular the case for the en- 
tropy dominated properties of membranes where 
assumptions of small configurational fluctuations 
or perturbative considerations fail. But even for a 
description of the membrane shapes at the mean- 
field level there are many challenges. An alterna- 
tive to the analytical approach is computer sim- 
ulations of self- avoiding fluid surfaces, which is 
viable both for studies of non- perturbative phe- 
nomena and shape transformations. The numeri- 
cal models of fluid membranes have been analyzed 
extensively, in particular plaquette models, where 
the surface is constituted by the plaquettes of a 
three-dimensional (3D) lattice [2HI], or 0(n) lat- 
tice gauge models for n — >• in 3D [5 . A drawback 
with the regular lattice based models of fluid mem- 
branes is the discrete nature of the surface config- 
urations, which make a detailed description of sur- 
face properties impossible and introduce phenom- 
ena which are not relevant for fluid membranes, 
e.g., the roughening transition. 
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The third class of numerical models for membranes 
is the triangulated random surfaces, which were 
introduced in statistical mechanics in context of 
Euclidean string theory[6]-[9]. Combined with sim- 
ulation techniques for self-avoiding polymers, the 
triangulated random surfaces served as models for 
lipid membrane conformations [10 . The fluid na- 
ture of the membrane is represented by a planar, 
triangular lattice structure, which is allowed to 
change connectivity throughout the simulation. A 
major advantage of these dynamically triangulated 
surface models is that discrete surface operators 
can be established which posses a simple contin- 
uum limit. The results from computer simulations 
of randomly triangulated surfaces can thus be in- 
terpreted in terms of continuum theory of mem- 
branes, the related literature has been reviewed in 

miiia. 

So far triangulated surface models only allowed 
for computer simulations of membranes equipped 
with pseudo scalar or scalar order parameters, e.g., 
mean curvature and density, while many interest- 
ing physical questions arises when vector or tensor 
order parameter fields are present in the plane of 
the membrane. For instance, tilting of the lipid 
molecules with respect to the surface normal, oc- 
curring in several of the ordered phases of lipid 
bilayers, give rise to in-plane orientational order- 
ing [13] . Furthermore, two good experimental ev- 
idences for the hexatic nature of the gel phase 
of lipid bilayer membranes have been reported 
recently [T4| H5]. Several classes of membrane in- 
clusions have the character of in-plane nemato- 
gens, e.g., antimicrobial peptides [16] and Bar do- 
main proteins, also see [17] and references within. 



In-plane orientational order in membranes has re- 
ceived major attention in the theoretical literature. 
In particular the properties of hexatic membranes 
[18] [19] and the Kosterlitz-Thouless transition phe- 
nomena on membranes [20 , the effect of lipid tilt 
and chirality [21-26 , and the effect of surfactant 
polar head order [27] . 

Here we present an approach to triangulated 
surface models of fluid membranes by combining 
the existing simulation technique of dynamical tri- 
angulation with an approach to compute the dis- 
cretized local curvature tensor. The properties of 
the random surface in the new description are con- 
sistent with those from earlier models. 

Furthermore, we study membranes with in-plane 
nematic order and show that it can give rise to non- 
trivial shapes. The paper is organized as follows: 
Sec. II introduces continuum models of mem- 
branes, the Helfrich Hamiltonian and its extension 
to include in-plane nematic fields with explicit cou- 
pling to the membrane curvature. In Sec. Ill we 
present the triangulated surface model which in- 
cludes a detailed description of the local surface to- 
pography, parallel transport along the surface and 
our numerical implementation of the model. The 
Monte Carlo procedure for computer simulation of 
the equilibrium properties of the triangulated sur- 
face model with in-plane orientational fields is de- 
scribed in Sec. IV. In Sec. V we characterize the 
nature of the triangulated surface for different val- 
ues of the bending moduli, without any in-plane 
order, and compare our results with that obtained 
from earlier simulations of membranes. In Sec. VI 
we discuss some examples, in our discretized mem- 
brane model, where the effects of the in-plane or- 
dering lead to some interesting shapes. 



II. CONTINUUM MODELS 

It has for long time been recognized that the 
large scale conformations of a simple closed fluid 
lipid membrane can be modeled by the Helfrich 
curvature energy functional [I] 

H c = | JdA (2M - 2C ) 2 + | JdA K (1) 

It is a purely geometrical model, where the char- 
acteristics of the surface is described by the con- 
formation of the membrane governed by the mate- 
rial constants, n the elastic bending rigidity, R the 
Gauss curvature modulus and Co the spontaneous 
mean curvature. K and M are the local Gauss and 
mean curvature of the surface respectively. There 
are several possible extensions of Eq.Q, e.g., de- 
scribing the effects of membrane inclusions , in- 
plane density fluctuations or in-plane order. Here 



we will discuss simple extensions of Eq.Q, now 
involving in-plane vector n or a nematic tensor or- 
dering field \ {h<&n). For a vector field, represented 
by an unit vector n, there is only one possible rele- 
vant extension of Eq.Q, to the lowest order in the 
order parameter [19], 
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which facilitates an implicit coupling of the mem- 
brane geometry to the ordering field. Ka is the 
stiffness constant and V is the covariant gradient. 
The model and its extensions have been analyzed 
in great detail (for review see chapters by Nelson, 
David and by Gompper and Kroll in [E]). For a 
nematic field, to the same order, the correspond- 
ing term is the well known Frank's free energy for 
nematics [281 
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h 1 - is orthogonal to h in the same plane. The in- 
plane Div(n) and Div(n^) are the splay and bend 
contributions of the nematic field, and K\ and K3 
are the corresponding Frank constants. The in- 
plane nematic field gives rise to a number of new 
relevant couplings between the ordering field and 
the curvature tensor [29]. A natural form of the 
free energy, that describes an explicit coupling be- 
tween the orientational field and the curvature ten- 
sor, is given by [2TH25| [30] 
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where, H n \\ is the directional curvature along ft 
and H n ± is the directional curvature along rr 1 . 

c|J and Cq are the corresponding spontaneous cur- 
vatures, ft || and k± respectively are the bending 
stiffness along h and h 1 - . 

III. TRIANGULATED SURFACE MODEL 

In this section we will consider discretized 
surfaces with the topology of a sphere, while the 
considerations can readily be extended to closed 
triangulated surfaces of arbitrary topology [3T| [32]. 
Contrary to the standard differential geometry of 
continuum models, the discretized formulation in 
this section is given in Cartesian coordinates. The 
surface is discretized by a triangulation T N con- 
sisting of N vertices connected by Nl = 3 (TV — 2) 
links, or tethers, forming closed planar graphs. 
The graph form a system of = 2 (TV — 2) 
triangles corresponding to a surface with total 
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FIG. 1: (Color online) Surface patch in a one ring 
neighborhood around vertex v. The edge e 
connects, in this description, v to 1. The edge 
vector is R{e) = X(l) - X(v) and N(e) is its 
normal. Edge e is shared by two faces fi(e) and 
/ 2 (e) with iV(/i(e)) and N(f 2 (e)), respectively, 
being their normals. The normal to vertex v is 
represented by N(v). 



Euler index x = X — N T — Nl = 2. Each vertex 
v takes a position X(v) in R 3 . The triangulation 
and the vertex position together form a discretized 
surface, a patch of which is given in Fig. [1] The 
self-avoidance of the surface is ensured by assign- 
ing a hard core spherical bead of unit diameter to 
each vertex and a maximal tether distance of y/3. 
This is in general not sufficient to impose strict 
self avoidance [33], [34], but a mild constraint on 
the dihedral angle between two faces sharing a 
tether restores self avoidance. 

The in-plane orientational field can be included 
by defining a unit vector n(v) in the tangent plane 
at each vertex v. In the following we will give 
meaning to this statement by analysis of the local 
surface topography and in turn calculate the cur- 
vature tensor, principal directions and curvature 
invariants [35] [36]. The approach is based on the 
construction of the discretized "shape operator" 
given by the differential form — dN in the plane of 
the surface, which contains all information about 
the local surface topography. 

Consider a local neighborhood around a vertex v, 
as shown in Fig. [l] R(e) is the edge vector that 
links v to a neighboring vertex. The set of edges 
linked to v is {e} vi while the oriented triangles or 
faces with v as one of their vertex is {f} v . The 
calculation of the surface quantifiers at v is re- 
stricted to the one ring neighborhood around it, 
which is well defined by {e} v and {f} v . Simi- 
larly the set of faces sharing an edge is given by 



{f}e — [/i( e )> /2(e)]. The normal to an edge e 
then is defined as, 



N(e) = 



N[fi(e)]+N[f 2 (e)] 



N[fi(e)}+N[f 2 (e)} 



(5) 



where N[fi(e)] and N[f2(e)] are the unit normal 
vectors to faces fi(e) and /2(e) respectively. 

We will now construct the shape operator at ev- 
ery vertex v. Toward this, we define 



H(e) = 2 R(e) 



cos 



*(e) 



(6) 



which quantifies the curvature contribution along 
the direction mutually perpendicular to R(e) and 
N(e) [35H37] . 3>(e) is the signed dihedral angle 
between the faces, j\ (e), /2(e), sharing the edge e 
calculated as 



$(e) = sign [{N[h(e)} x N[f 2 (e)}} ■ R(e) 
[N[fi(e)]-N[f 2 (e)]\ + tt. 



(7) 
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The discretized "shape operator" , which quantifies 
both the curvature and the orientation of e is thus 
the tensor 

Se(e) = H(e) \R(e) x N(e)\ \R(e) x N(e) \ , (8) 

where R(e) = R(e)/\R(e)\ is the unit vector 
along edge e. Having defined the shape operators, 
{Se(e)}, along the edges of the vertex v, we now 
proceed to compute the shape operator at v. The 
normal to the surface at v can be calculated as, 
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E if}v m(f)]N(f) 



(9) 



with A(f) denoting the surface area of the face / 
and the normalized weight factor fi[y4(/)] is pro- 
portional to the area of the face. The projection 
operator, P(v) = 1 — N(v)N(v), projects {S e (e)} 
on to the tangent plane at v [35, 36 . The shape 
operator at the vertex v is then a weighted sum of 
these projections given by 

§vW = 7L E^W^^EW- ( 10 ) 

1 ^ {eh 

A(v) — ^Z{f} v is the average surface area 

around v, while the weight factor for an edge is 
calculated as W(e) = N(v) • N(e). The shape op- 
erator Eq.([lO|) at the vertex v is expressed in co- 
ordinates of the global reference system [x,y,z]. 
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Notice that, by construction, the vertex normal 
N(v) is an eigenvector of, S v (w), corresponding 
to eigenvalue zero. The two other principal direc- 
tions ti(v), i 2 (v), whose corresponding eigenval- 
ues are the principal curvatures, define the tangent 
plane at the vertex v. A local coordinate frame 
called the Darboux frame [ii(v), N(v)], see 

Fig(2j can then be defined at v. The transfor- 
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FIG. 2: (Color online) Transformation from a 
global to local coordinate frame. 

mation from the global to Darboux frame, see 
Fig. [2j is obtained by first applying a Householder 
transformation(H), see Appendix A, to rotate the 
global z direction into N(y), while x and y are 
rotated into vectors x' \y' in the tangent plane at 
the vertex v. The shape operator, at v, in this 
frame C(v) = Hl(v) S v (w) H(v) is a 2x2 minor, 
with the two principal curvatures c\{v) and C2(v) 
as its eigenvalues. The corresponding eigenvector 
matrix E(v) transform [x , y , N(v)] into the Dar- 
boux frame at v. Any vector in the global frame, 
can now be transformed to this local frame by the 
transformation matrix E(v)H(i;). 

We are now in the position to write up the dis- 
cretized form of Helfrich's free energy, at a vertex 
v, based on the local curvature invariant M(v) = 
[ci(v) + c 2 (v)}/2 and K(v) = 2a(v)c 2 (v): 
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(11) 

The calculation of the discrete curvature tensor has 
been performed by other methods [38[ [39], how- 
ever we find the method used in this paper to be 
the most accurate in describing surfaces with pre- 
scribed geometry. The local Darboux frame is very 
useful for the characterization of an in-plane vec- 
tor field h(v). For convenience, we choose c\(v) to 
be the maximum principal curvature and i\(v) the 
corresponding principal direction. The local orien- 
tational angle cp(v) of n(y) will always refer to this 
Darboux frame. 

In order to compare the orientation of two distant 
in-plane vectors at the surface, it is necessary to 
perform parallel transport of the vectors on the dis- 
cretized surface. In practice we need only to define 
the parallel transport between neighboring ver- 
tices, i.e. a transformation h(v ) —> T(v, v )n(v), 



which brings h(v) correctly into the tangent plane 
of the vertex v , so that its angle with respect to 
the geodesic connecting v and v is preserved. If 
f(v,v ) is the unit vector connecting a vertex v to 
its neighbor v and ((v)=~P(v)r(v,v ) and ((v ) = 
P_(v )r(v ,v) are its projection on to the tangent 
planes at v and v ; then our best estimate for 
the directions of the geodesic connecting them, are 
the unit vectors C(v), ((v ). The decomposition of 
h(v) along the orientation of the geodesic and its 
perpendicular in the tangent plane of v is thus: 



n(y) = h(v) ■ £(v) C(v)+ 
h(v) • (N(v) x C(v)) (n(v) x 



(12) 



Parallelism now demand that these coordinates, 
with respect to the geodesic orientation, are the 
same in the tangent plane of v , therefore: 



T(v,v)h(v) = [n(v) -C(v) ((v')+ 
{h(v)-(N(v)x((v))} [n(v')xC(v) 



(13) 



This parallel transport operation allow us to define 
the angle v ) between vectors in the tangent 
plane at neighboring vertices, and in turn their co- 
sine and sine as: 



cos((/)(v,v )) = ii(v ) • T(v, v )h(v); 



(14) 



sm((/)(v,v)) = N(v ) x h(v ) -T(v,v)n(v) 

We can now define the lattice models, correspond- 
ing to Eqs.([2| and ([3|, for the in-plane orienta- 
tional field, e.g., the XY- model on a random sur- 
face: 



n XY = 
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(15) 



or the Lebwohl-Lasher model on a random surface: 
cll f3 



LL 



2 E {^os'i^v'))-^ (16) 

(vv ) 



Furthermore, we are now in a position to calculate, 
at a given vertex v, the directional curvatures along 
and perpendicular to the orientation of the in plane 
vector field h(v) by use of Gauss formula: 

M(v) || = ci (v) cos 2 [(p(v)] + c 2 (v) sin 2 [<p(v)] 
M(v)± = ci(v) sin 2 [(/?(v)] + c 2 (v) cos 2 [(/p('u)] 

(17) 
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IV. MONTE CARLO PROCEDURE 

The equilibrium properties of the discretized 
surface can now be evaluated from the analysis of 
the total partition function, i.e., the sum of Boltz- 
mann factors for all surface configurations and tri- 
angulations. For simplicity, we consider the situa- 
tion with just one in-plane orientational h(v) field 
defined at each vertex 

N 

Z(N, k, R, e, ..) = ^ J2 II / djl ^ I Mv) 

exp (h c ({X},T n ,W}) + Usas)) (18) 

where Usas is the potential that ensures the 
self- avoidance of the surface and cp(v) is integrated 
over the unit circle or half unit circle for the XY 
field and the nematic field respectively. {X} and 
{ip} are, respectively, the complete set of vertex 
positions and orientational angles. Further, we set 
P = j^pf = 1. In practice, a surface configuration 

is represented by a tuple r] = ({X}, T N , {</?}), 
which must be updated during the Monte Carlo 
simulation procedure. The Monte Carlo updating 
scheme can now be decomposed into three move 
classes, so each of the three sets of degrees of 
freedom are updated independently to keep it 
simple and ensure fulfillment of detailed balance: 

Vertex shifts: represent the updates of the 
vertex positions, keeping T N ,{(p} fixed, thus 
allowing for shape changes of the membrane. 
The attempt probability to change to a new 
configuration rj' = ({A'}, T N , {<£}), with a chosen 
vertex moved to a new position within a cube 
of side 2<j centered around its old position, is 
u(rj\rj') = u(r}'\rj) = ([ 2cr ] 3 ^) -1 - & is appro- 
priately chosen to get a reasonable acceptance 
rate of 30-50%. In our simulations a = 0.1 is 
chosen. With this surface updating operation, 
the curvature tensor and thus the principal axis 
changes. Since the angle {(/?} is kept fixed, the 
set of orientations {n}, in the global frame, are 
changed following the local surface configuration, 
Fig. ga). 

Link flip: represents updating of the trian- 
gulation. Here a link, e connecting a vertex v 
to v , is picked at random and an attempt is 
made to flip it to the pair of opposite vertices 
common to v and v . The attempt probability 
to change to configuration rj' = ({-?}, T /Ar , {(/?}) 
is then Lj(r)\rf) = uj(r]'\r]) = 1/Nl- Similar to 
vertex shifts, the actual orientations {h} are now 
changed, following the local surface configuration, 
Fig. gb). 

Angle rotation: the orientation of the in-plane 




FIG. 3: (Color online) Monte Carlo moves, a) 
vertex shift, b) link flip and c) angle rotation. 
Surface vector field, n, is represented with solid 
arrow while principal directions e\ and &2 are 
marked with dotted arrows, ip is the angle n 
subtends with e\. 



vector h(v), at a randomly chosen vertex v, 
is updated. The vector is rotated to a new, 
randomly chosen, direction in the tangent plane, 
keeping the vertex positions and link directions 
fixed. As a result of which the orientational 
angle is now ip (v) = <p(v) + Acp(v). The attempt 
probability to configuration v[ — ({X},T N ,{<p }) 
is u(r]\r] ) = (jj(r}'\rj) = (2cr^7V) _1 , where <C tt is 
the maximum increment of the angle. The surface 
topography is not affected by this move, Fig. [3^c). 

For each of the above moves, the acceptance 
probability is: 

accfrW) = min(l, ^-M exp(-/3 (h( V ') - H(rj))) 
uj{r]\r] ) V / 

(19) 

The duration of a Monte Carlo simulation is 

measured in MCS (Monte Carlo sweeps per Site), 
which represents N attempted vertex moves, 3(N— 
2) attempted flips and N attempted rotations of ft. 
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V. RESULTS AND DISCUSSION 

A. Vesicles with no in-plane order 

In the first part of this section we will discuss the 
properties of this new discretized random surface 
description of membrane conformations for a sim- 
ple, closed, fluid membrane of spherical topology, 
with no in-plane order. All simulations reported in 
this paper are carried out with ft = 0. System sizes 
in the range N = 77 to 3677 and bending rigidity 
in the range ft = to 1000 were investigated to 
compare it with the previously known results for 
these systems. Applying the equipartition theo- 
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FIG. 4: (Color online) versus ft for stiff 
membranes and varying N . Note that for large 
values of ft the line approaches 47r, the value of 
(H c ) over bending modulus, for a smooth sphere 
with bending modulus ft. Inset shows a time 
series of the curvature energy for system of 
N = 2030 vertices for k = 200. 



limit ft* and C max , in dimensionless units, sat- 
urates to approximately 4.4 and 1.4 respectively. 
The smooth and finite nature of C(N, ft) for large 
N shows that this measure does not indicate the 
presence of a first order or a continuous transition 
in the thermodynamic limit. However, a continu- 
ous transition cannot be completely ruled out. If ft 
is an irrelevant thermodynamic variable under RG 
transformation, it just leaves a cusp in C(iV, ft) at 
the transition, a similar phenomena is well-known 
for the A-transition of He 3 -He 4 mixtures [43 . Note 
that the value of ft* appears to be roughly five 
times that of the previously reported values [40) [42]. 
This is a clear indication of that the new measure 
of local mean curvature differs from that used pre- 
viously, although the prediction of a low ft cusp 
in C(7V, ft) persists. A simple quantifier of mem- 




FIG. 5: (Color online) Specific heat C(N, ft) 

versus ft for varying N. The position of 
C max (7V, n)(circles fitted with solid line) and 
ft* (N) (squares fitted with dotted line) are shown 
in the inset. 



rem to Gaussian or quasi-spherical configurational 
fluctuations shows that the expected behavior is 
^ — > 8tt + ^f 1 ^. In Fig. [4J it is shown that 
the ensemble averaged curvature energy of the vesi- 
cle, indeed approaches Sir for large ft. In the 
opposite limit of small ft the literature is largely 
focused on the crumpling transition. Such a tran- 
sition should be indicated by the presence of a peak 
or a cusp in the specific heat, 

C(N,k) = ±((H 2 c )-(H c ) 2 ). (20) 

C(N , ft) , as a function of ft for different TV, is shown 
in Fig. [5] The shape of the curve is similar to what 
has been reported by earlier simulations [40ti42] , 
As reported in these papers, we find that the peak 
height ((7 max ) stops growing and the peak position 
(ft*) saturates to a constant value beyond system 
size N « 500 ( see Fig. [5]). In the aysmptotic 



brane conformations used in triangulated surface 
simulations is the gyration tensor 

1 N 

v,v 

with Rq = Tr(G) as the simplest invariant. For 
the flexible, tethered, self-avoiding random sur- 
faces Rq ~ N a . Earlier simulations report a = 
0.8 [34] and a = 1. [40 . As shown in Fig. [6] we find 
that Rq ex N for all values of ft, which is charac- 
teristic of the self-avoiding branched polymer and 
quasi spherical configurations [4] . The similarity of 
the exponent makes an analysis of the cross-over, 
between the branched polymer configurations at 
low ft and quasi- spherical shapes at high ft, very 
difficult by use of Rq. This is better accomplished 
by analysis of the vesicle volume, which in previ- 
ous vesicle simulations have been shown to obey 
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FIG. 6: (Color online) linear scaling of R g as a 
function of N for various n is shown. Entropic 
domination in lower k regime brings in a large 
spread in the values of R 2 g /N for k < 1.0 and are 
not shown here. 



the simple scaling ansatz V = AT 2 f ^ aN /£ p (k) 

where f(x) is a scaling function and £ p (k) is a 
cross-over length scale, identified with the persis- 
tence length [33j [44]. This universal scaling be- 
havior also holds for our new triangulated surface 
model as shown in Fig. [7| Here, for each £ p is 



plots shown in Fig. [7j predicts a different depen- 
dence on n ( see Fig. [8|. In the flexible regime, 
k < 3/C5T, an approximate exponential behavior 
exp(c/c/fcbT), c ~ 7r/6 is seen, while in the semi- 
flexible regime, n > a stronger dependence 
of £, p on k, is found. Our data does not allow for a 
determination of the asymptotic behavior of £ p (k) 
for large k. The scaling function where 
\I/ = x/iV^" 1 , is a constant for small \I/ ( semi- 
flexible regime ) and is ~ for large values of 
^, indicating a branched polymer behavior in the 
flexible regime, see Fig. [7[ This suggest that, in 
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FIG. 8: (Color online) Persistence length, ^, as a 
function of n 
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FIG. 7: (Color online) Universal scaling function 
describing the dependence of volume on the 
system size. The data collapse is obtained by 
determining £ p for each k separately. 

chosen such that we obtain good data collapse. It 
has been found by RG-analysis that £ p for a fluctu- 
ating smooth continuous surface, embedded in 3D 
space, depends on n as ex${4n k / SK^T) [45]. This 
dependence has been verified numerically by previ- 
ous triangulated surface simulations [33]. However, 
the persistence length, obtained from the scaling 



this model, the effective bending rigidity is a de- 
creasing function of temperature, with c saturat- 
ing to 47r/3 at low temperatures. Overall, we have 
shown in this section that the new algorithm repro- 
duce the expected behavior of vesicles governed by 
Helfrich's free energy , given in Eq.Q, in the rigid 
regime of high k. In the flexible to semi-flexible 
regimes of low n values, our new numerical repre- 
sentation of the geometry and energetics of vesicles 
produce a behavior which is qualitatively in agree- 
ment with previous triangulated surface models of 
vesicles. However, the cusp in the specific heat has 
shifted to higher n value. The flexible regime at 
k values below the cusp is more rigid compared to 
previous models with an approximate exponential 
dependence between the persistence length and n 
and i p (n*) — 10. Above the increase in £ p is 
much stronger. We attribute the differences be- 
tween the present model and previous models to 
the use of different surface quantifiers. 



B. Membranes with in-plane nematic order 

We will consider the case of a randomly tri- 
angulated surface with an in-plane nematic field. 
These systems, in the continuum limit, are de- 
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FIG. 9: (Color online) Equilibrium configuration of a nematic embedded vesicle with k = 0, cjj = 0, 
cll = 3.0, k±_ = and (a) k\\ = 0, mere presence of an nematic field in the ordered phase cuts off the 
entropy dominated branched polymer phase seen otherwise (b) «y =20 and (c) a corner with defect of 
index + | is shown for K\\ = 0. All data are for a triangulated surface with 1202 vertices. 



scribed by a free energy functional which contains, 
in addition to the basic Helfrich curvature elastic 
part Eq.Q, terms describing nematic-nematic in- 
teractions Eq.([3| and the coupling of the nematic 
field to the membrane curvature Eq.Q. For the 
discretized nematic-nematic interactions we have 
employed the Lebwohl-Lasher[46j H7] model, de- 
scribed in Eq.(16), which corresponds the one con- 



stant approximation of Frank's free energy given 
in Eq.(([3|). The total discretized free energy func- 
tional then takes the form 



V = l 



E E { 

v v'e{v} 

N 

yJ2 [ h hv),w - 



3 / 1 

-cos 2 (cj)(v,v )) - - 



A(v) 



[ H nM t ± -C^Y A{v), 



v=l 
N 



(21) 



where, #n(»,|| = n\{v) c\(v) + n 2 (v) c 2 (v) and 
Hh(v),± = ni(v) 2 C2(v) + ri2(v) 2 ci(v) are the di- 
rectional curvatures at a vertex v, see Eq.(17). 
M(v) = [c\(v)+C2(v)\/2 is the corresponding mean 
curvature. Note that this free energy is expressed 
in the local Darboux frame of reference, described 
in Sec. 



Ill 



774 (v) and ri2(v) are the components of 
the nematic director in this local frame, and c\(v) 
and C2(v) are its principal curvatures. A(v) is the 
area of the polygonal surface defined by its nearest 
neighbors. 

We will, in what follows, demonstrate the use of 
the algorithm by studying the effect of in-plane ori- 
ent ational ordering on membrane conformations. 



A detailed quantitative analysis and phase dia- 
gram of the vesicles shapes and in-plane ordering 
that can result from Eq.(21 ) will be published else- 
where. 



1. Membrane stiffness originating from the ne- 
matic field 

First we consider the case with k = 0, k\\ ^ 
and c\\ =0. We choose k± = so that the nematic 
field does not directly influence the bending modu- 
lus perpendicular to it. Such a situation may arise 
in the case of long thread like inclusions, cll = 3 
is chosen to favor in-plane nematic order. 

Characteristic equilibrium configurations corre- 
sponding to fty = and 20 are shown in Fig.[9j For 
ft || =0 the common shapes are deformed tetrahe- 
drons with four well-defined corner points. The in- 
plane orientational field displays perfect nematic 
ordering except at the corner points where a discli- 
nation with index 1/2 is located. A snapshot of one 
of these disclinations is shown in Fig. |9jc. Since 
these are the only disclinations, the total index is 2, 
in accordance with Poincare's index theorem. The 
surface appear crinkled with local scale roughness. 
For n || =20 the vesicle shape becomes elongated, 
with the long axis following the orientation of the 
nematic field, with sharp ends. The two 1/2 de- 
fects are now joined to form a defect of index 1, 
and is located at the two ends. 

Membrane without stiffness and nematic degrees 
of freedom has branched polymer configurations. 
While our simulations show that, for the same sys- 
tem size, such a phase is absent in membranes with 
in-plane order. It thus follows that in-plane or- 
dering induces configur ational stiffness of vesicles. 
Signature of this stiffness can also be seen in the 
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distribution of eigenvalues of the gyration tensor 
for two different values of €ll- As can be seen in 
Fig. [lOj the distribution of higher eigenvalues are 
narrower for higher cll, indicating stiffening. We 
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FIG. 10: (Color online) Distribution of the 
eigenvalues (A^) of the gyration tensor, such that 
Ai < A 2 < A 3 , for k = k\\= k± = 0, (a) e LL = 1 
and (b) cll = 10. 

note that the anisotropic elasticity of the mem- 
brane, arising through this nematic orientation, is 
similar to that suggested by Fosnaric et al [48] . 




FIG. 11: (Color online) Configurations of 
membranes, with k,\\ = 20, cjj = 0.5,6/,/, = 3, 
k ± = for (a) k = 2.5 and (b)tt = 10 



Positive spontaneous curvature 



Making c[| > imply that the nematic field fa- 
vors a specific value of positive curvature along the 



direction of its axis. In Fig. 11 is shown representa- 
tive equilibrium configurations for cll — 3, k± — 
0, k\\ = 20, cjj = 0.5 and k = 2.5(a), 10(b). For 
k = 2.5 the vesicle shape transforms to branched 
structure with long irregular tubes of varying ra- 
dius. The nematic field now spirals around the 
tubes. The angle made by the nematic field with 
the azimuthal direction increases with decrease in 
local tube radius. The caps of the tubular struc- 
tures are quipped with disclination pairs of index 
1/2, while two disclinations with index —1/2 are 
situated in the branchpoints. The tubes them- 
selves tend to spiral over longer distances, as can 
be seen from Fig. [TTJa. This spiraling can both be 
right and left handed, indicating no chiral prefer- 
ence. For k = 10 this picture persists, except the 
nematic ordering match up with the azimuthal di- 
rection of the tubes, no chiral ordering of the tubes 
are observed and the tube radius match with the 
that set by cjj. 



Negative spontaneous curvature 



Negative spontaneous curvature, Cq < 0, implies 
that the nematic will now prefer to orient along 
directions where the membrane curvature is nega- 
tive ( curved into the vesicle). In Fig. 



12 



is shown 

examples of equilibrium configurations for n = 
and k = 10, where cll = 3.0, fty = 30, k±_ = 0, 
cjj = —0.5 and n±_ = 0. Inward tubulation results 
in stiffening of the outer boundary of the vesicle 



as shown in Fig. 12 In contrary to the tubulation 



seen in the case of Cq > 0, self avoidance condi- 
tion of the membrane now prevents complete tube 
formation. Similar to the cjj > case, increasing 
k, increases the thickness of the tubes. The ne- 
matic ordering is along the azimuthal direction for 
k, = 10, while spiraling is observable for n = 0. On 
the outer surface, defects of index — \ are clearly 
observed. 
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FIG. 12: (Color online) Configurations of membranes, with fty = 30, cll = 3.0, cjj = —0.5, ft^ = for 
(a) ft = and (b) ft = 10. (c) is the mesh representation of the surface in (b) which clearly shows tubes 

grown inwards. 



VI. CONCLUSION 

We have presented a methodology for calcu- 
lating surface quantifiers on a self-avoiding tri- 
angulated random surface models of fluid mem- 
branes. The method involves calculations of the 
local geometrical properties at the vertex positions 
of the surface, e.g., calculation of the local Dar- 
boux frame and the principal curvature radii of 
the surface. We have described a procedure for 
parallel transport of in-plane vectors between ver- 
tex points. We have implemented the numerical 
model and performed Monte Carlo simulations of 
the equilibrium properties of the surface. The sim- 
ulations of the discretized form for the Helfrich's 
free energy are in good qualitative agreement with 
the results from previous numerical simulations. In 
the flexible limit of low bending rigidity the mem- 
brane scales as a branched polymer and a scaling 
relation involving volume, system size and persis- 
tence length holds. For small values of ft, calcu- 
lations using the new discrete Hamiltonian shows 
a faster increase, as a function of ft, in the persis- 
tence length compared to the previous model. 

The model has been extended to include an in- 
plane nematic field and equilibrium shapes have 
been obtained for some simple examples. We show 
that the presence of a nematic ordering leads to 
suppression of the branched polymer phase even 
when the bare bending rigidity is zero. The con- 
formational changes in a fluid membrane brought 
about by the anisotropy in the bending rigidity and 
the spontaneous curvature induced by the nematic 
field are demonstrated. We have demonstrated 
that the presence of the in-plane nematic field leads 



to coupling between geometry and nematic defect 
structures of the membrane. It is shown that this 
coupling can lead to chiral structures in membrane 
even in the absence of explicit chiral terms in the 
Hamiltonian. 



Acknowledgements 

The MEMPHYS-Center for Biomembrane 
Physics is supported by the Danish National 
Research Foundation. Computations were carried 
out at the HPC facility at IIT Madras and Danish 
Center of Scientific Computing at SDU. 

Appendix A: Householder Transformation 

Consider two orthonormal frames of reference 
given by the coordinates (x , z) and (a, 6, c). The 
Householder matrix H, can be used to rotate z 
in frame 1 to c in frame 2, such that (£, y) now 
are some arbitrary vectors in the plane formed by 
(a, b). Define a vector, 



with a minus sign if 1 1 x — c\ | > | \x + c\ | and a plus if 
otherwise. The Householder matrix is then defined 

as, 

H = 1 - 2WW ] (A2) 
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